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ABSTRACT 

We perform a combined fit to angular power spectra of unresolved infrared (IR) point sources from the Planck 
satellite (at 217, 353, 545, and 857 GHz, over angular scales 100 < l < 2200), the Balloon-borne Large-Aperture 
Submillimeter Telescope (BLAST; 250, 350, and 500 pm\ 1000 < i < 9000), and from correlating BLAST 
and Atacama Cosmology Telescope (ACT; 148 and 218 GHz) maps. We find that the clustered power over the 
range of angular scales and frequencies considered is well fitted by a simple power law of the form C f clusl oc t~ n 
with n = 1.25 ± 0.06. While the IR sources are understood to lie at a range of redshifts, with a variety of 
dust properties, we find that the frequency dependence of the clustering power can be described by the square of 
a modified blackbody, B(v, T e g), with a single emissivity index ft — 2.20 ± 0.07 and effective temperature 
T e ff = 9.7 K. Our predictions for the clustering amplitude are consistent with existing ACT and South Pole 
Telescope results at around 150 and 220 GHz, as is our prediction for the effective dust spectral index, which we 
find to be a \ 50-220 = 3.68 ± 0.07 between 150 and 220 GHz. Our constraints on the clustering shape and frequency 
dependence can be used to model the IR clustering as a contaminant in cosmic microwave background anisotropy 
measurements. The combined Planck and BLAST data also rule out a linear bias clustering model. 

Key words: cosmic background radiation - cosmology: observations - infrared: diffuse background - infrared: 
galaxies - submillimeter: diffuse background 
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1. INTRODUCTION 

The angular power spectrum of cosmic microwave back- 

ground (CMB) temperature fluctuations currently provides 

vital constraints on cosmological models (e.g., Komatsu et al. 

2011; Larson et al. 2011). Experiments including the Atacama 

Cosmology Telescope (ACT), South Pole Telescope (SPT), and 
Planck satellite are now probing the CMB temperature power 

spectrum on arcminute scales (Das et al. 2011; Keisler et al. 

2011; Planck Collaboration et al. 20 1 1 a) . An improved measure- 

ment of the Silk damping tail (Silk 1968) improves constraints 

on, for instance, the scale dependence of primordial fluctua- 
tions, important for testing inflationary models, the number of 
relativistic species, and early-universe exotica (e.g., Komatsu 
et al. 2011). On arcminute scales the contribution to the angular 
power spectrum from the primary CMB fluctuations becomes 

subdominant to extragalactic foregrounds including infrared and 
radio point sources (e.g., White & Majumdar 2004; Righi et al. 
2008), and the thermal and kinetic Sunyaev-Zel’dovich effects 
(SZ; Sunyaev & Zel’dovich 1970). Extracting the CMB signal 
requires understanding how the contribution from these compo- 
nents varies with frequency and angular scale. 

The infrared (IR) point-source foreground is understood to 
originate from high-redshift (z ~ 1 — 4) star-forming galaxies 


whose rest-frame emission peaks in the far-infrared due to 
thermal emission from dust grains illuminated by starlight (e.g., 
Bond et al. 1986, 1991; Hughes et al. 1998; Blain et al. 1999; 
Draine 2003). While our work concerns observations made in 
the millimeter and submillimeter we refer to these dusty sources 
throughout as “IR” sources. Thermal dust emission from star- 
forming galaxies is also an important component of the cosmic 
infrared background (CIB — e.g., Puget et al. 1996). Galaxies 
trace the large-scale structure and so are clustered, with a scale- 
dependent contribution to the power spectrum in addition to 
Poisson shot noise (Peebles 1980). Clustering of IR sources 
was therefore expected (e.g., Bond 1996, and references therein; 
Negrello et al. 2007) and has been detected in Spitzer Space 
Telescope data at 160 pm (Lagache et al. 2007), by the Balloon- 
borne Large-Aperture Submillimeter Telescope (BLAST) and 
Herschel Space Obser\>atory at 250, 350, and 500 pm (Viero 
et al. 2009; Cooray et al. 2010; Amblard et al. 2011), in 
the microwave sky at around 150 and 220 GHz by SPT and 
ACT (Hall et al. 2010; Dunkley et al. 2011; Shirokoff et al. 
2011), and in early data from Planck (Planck Collaboration 
et al. 2011b, hereafter PI 1 ) . Correlations between clustering 
at different frequencies have also been detected: Hajian et al. 
(2012, hereafter HI 2) measure significant levels of correlation 
between BLAST maps and ACT maps at 148 and 218 GHz, 
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detecting a clustered component at 4a . PI 1 also find a significant 
correlation between the Planck 217 GHz maps and those at 
higher frequencies (353, 545, and 857 GHz). Studying the 
clustered power of dusty galaxies in the submillimeter is simpler 
than at the millimeter CMB bands because they are the dominant 
extragalactic signal. Knowing that significant overlap exists 
between the submillimeter galaxy population and those sources 
responsible for the IR foreground in the CMB bands suggests 
that information about the former can help our understanding of 
the latter. 

In this work, we combine large- and small-scale power spectra 
from Planck, BLAST, and correlations between BLAST and 
ACT to estimate the amplitude and scale dependence of the 
angular power spectrum of clustered IR sources. We present 
a simple power-law template to model the clustered source 
contribution that may be marginalized over when estimating 
cosmological parameters from foreground-contaminated CMB 
maps as in, for instance, Hall et al. (2010), Dunkley et al. (201 1), 
and Keisler et al. (201 1). 

Previous ACT and SPT results (Dunkley et al. 2011; 
Shirokoff et al. 2011) have found that the parameters extracted 
from their CMB spectra are not particularly dependent on the 
model adopted for the IR clustered power (see also Sehgal et al. 
2010; Fowler et al. 2010). The ACT and SPT data sets are not yet 
complete; the final data will include more sky coverage as well 
as measurements from additional frequency channels. Millea 
et al. (2012) find that modeling the IR point-source clustering 
incorrectly for the final combined Planck and ACT/SPT data 
sets could introduce a significant bias in cosmological param- 
eters (they estimate 1 a based on the discrepancy between two 
different IR clustering models). It is therefore important to un- 
derstand the scale and frequency dependence of the clustered 
power in preparation for this future analysis. Improving con- 
straints on IR clustering will also help constrain the SZ power 
spectrum. 

In this work, we are primarily concerned with the IR sources 
as a contaminant in CMB maps. While the physical properties 
of this high-redshift star-forming population are important for 
understanding the star formation history and galaxy evolution, 
in this analysis we do not attempt to extract information about, 
for example, the redshift distribution of the sources or the dark 
matter halos they occupy. 

In Section 2 we describe the data we use for our fitting, in 
Section 3 we explain our assumptions and methods, results are 
presented in Section 4, and a conclusion follows in Section 5. 

2. DATA 

Throughout this work we use “auto-spectrum” to refer to 
a power spectrum calculated by correlating two maps at the 
same frequency, and “cross-spectrum” to refer to a spec- 
trum calculated by cross-correlating maps at different frequen- 
cies. “BLAST x ACT” is to be understood to refer to the 
BLAST/ACT cross-spectra and so on. 

We use the BLAST 250, 350, and 500 jim auto-spectra, 250 x 
350, 250 x 500, and 350 x 500 /im cross-spectra, and BLAST 
250, 350, and 500 /zm x ACT 148 and 218 GHz cross-spectra 
from H12, and Planck 857, 545, 353, and 217 GHz auto-spectra 
from Pll to construct the template. The ACT data are from 
the 2008 observing season and the BLAST x ACT spectra 
were calculated from the ~8.6 deg 2 common to both sets of 
maps, as described in H12. We take the quadrature sum of the 
statistical and beam systematic uncertainties in Table 4 of PI 1 
as the error on each Planck data point, neglecting any possible 


correlation in the beam uncertainty across different angular 
scales. We find that allowing such a correlation has minimal 
effect on our results (Section 4.1.3). The BLAST and ACT 
beam uncertainties are subdominant to the statistical (noise) 
uncertainty across the range of angular scales covered by the 
BLAST and BLAST x ACT data and so we likewise neglect any 
correlation in the errors on these spectra. The temperature-to- 
flux conversion factors given in Table 4 of PI 1 assume a source 
spectral energy distribution (SED) that varies as /(v ) oc v -1 ; a 
correction is then applied to convert to the real flux units (see 
Section 5.5 of PI 1 and Section 7.4.2 of Planck HFI Core Team 
et al. 2011). 

We subtract the Galactic dust emission (cirrus) component 
in the BLAST and BLAST x ACT spectra as described in 
Section 4.3 of HI 2, assuming the cirrus contribution varies with 
angular scale as i~ 2n . These spectra are not very sensitive to the 
details of the cirrus treatment, since they are from a relatively 
cirrus-free patch of sky (see also Das et al. 201 1). We assume 
that any contribution to the power spectra other than that from IR 
galaxies (for instance, radio-galaxy-IR or SZ-IR correlations) 
is negligible. Removal of power in the Planck spectra that is not 
from extragalactic IR point sources is described in Sections 2 
and 3 of Pll. 

The Planck and BLAST auto-spectra are presented in 
Figure 1. We also show the 218 GHz ACT power spectrum from 
Das et al. (201 1) and the 220 GHz SPT spectrum from Shirokoff 
et al. (201 1). A WMAP-1 best-fit ACDM CMB power spectrum 
(Komatsu et al. 2011) has been subtracted from these points; 
we have not subtracted any power from radio point sources or 
the kinetic SZ effect, as these contributions are likely to be sub- 
dominant (Dunkley et al. 2011; Shirokoff et al. 2011). All the 
spectra show similar angular scale dependence despite spanning 
a broad range of frequencies and almost two decades in angular 
scale. Note that a color correction to account for the different 
bandpass filters is required to make the Planck and BLAST 
spectra directly comparable — see Section 4 for details. 

The Planck, BLAST, and ACT maps are calibrated using com- 
parisons to various measurements including the orbital CMB 
dipole, FIRAS data, and planetary temperature. Uncertainty in 
these calibrations must be accounted for in order to perform 
joint fits to the power spectra extracted from the different maps. 
The absolute photometric calibration uncertainties are 7% for 
the Planck 857 and 545 GHz maps, 2% for the Planck 353 and 
217 GHz maps (see Pll and Planck HFI Core Team et al. 2011), 
9.5%, 8.7%, and 9.2% for the BLAST 250, 350, and 500 /im 
maps (Truch et al. 2009), 7% for the ACT 218 GHz map, and 
2% for the ACT 148 GHz map (Hajian et al. 2011). Further- 
more, the uncertainties in the BLAST calibrations are highly 
correlated; this is because the dominant source of uncertainty is 
the SED of the star used for the calibration of all three bands 
(Truch et al. 2009). 

3. POWER-LAW CLUSTERING TEMPLATE 

In the halo model formalism (e.g.. Peacock & Smith 2000; 
Scoccimarro et al. 2001; Cooray & Sheth 2002; see also 
earlier work, e.g.. Bond 1996, and references therein), the 
two-point function of galaxy clustering is dominated on large 
scales by pairs of sources in different dark matter halos (the 
“two-halo” term — roughly corresponding to the linear clus- 
tering regime), while on small scales galaxy pairs occupying 
the same halo (“one-halo” term — nonlinear regime) become 
dominant. Although these two components have different 
scale dependence, an angular correlation function varying as 
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Figure 1. Planck (857, 545, 353, and 217 GHz) and BLAST (250, 350, and 500 fim) IR point-source power spectra. The spectra include both shot noise and clustered 
components. While these data span a broad range of frequency and angular scale there is a notable similarity in the angular scale dependence. The Planck error bars 
include both the statistical uncertainties and estimates of the systematic beam uncertainty given in Table 4 of Planck Collaboration et al. (2011b). Data points at 
i > 2000 from ACT at 218 GHz (Das et al. 2011) and SPT at 220 GHz (Shirokoff et al. 2011) have been included for comparison with the Planck 217 GHz spectrum. 
We have subtracted a WMAP-7 best-fit ACDM CMB component from these spectra. No corrections for the different bandpass filter profiles have been applied; when 
the different filters and photometric calibration uncertainties are accounted for the Planck and BLAST data are in good agreement (Section 4). 

(A color version of this figure is available in the online journal.) 


w(0) oc 9~ s with S ~ 0.8 (Peebles 1980, corresponding to 
clustering power C| lust oc A 1 ' 2 ) has been found to adequately 
describe the clustering of, for instance, Lyman break galaxies at 
z ~ 3 (Giavalisco et al. 1998) and local SDSS galaxies (Zehavi 
et al. 2002). Significant deviations from power-law behavior 
have however been observed in analyses of more recent and 
comprehensive SDSS data (e.g., Zehavi et al. 2004, 2005; Blake 
et al. 2008), as well as in high-redshift galaxies (Ouchi et al. 
2005; Lee et al. 2006; Coil et al. 2006; Wake et al. 2011). For 
a recent discussion of the physical origins of power-law galaxy 
correlation functions see Watson et al. (2011). 

PI 1 and Amblard et al. (201 1) find that a power law provides 
an adequate fit to the existing Planck and Herschel /SPIRE 
unresolved IR point-source spectra (see Tables 6 and S 1 in those 
papers, respectively). This would suggest that IR spectra are not 
yet of sufficient quality to reveal deviations from power-law 
clustering behavior and we therefore also adopt a power law for 
the clustering component in this work. 

We model the total IR point-source power spectrum from 
correlating maps at frequencies v\ and v 2 (vi = v 2 for the 
auto-spectra) as 

/ l 

Ci(v u v 2 ) = A c (v u v 2 )l— J +C P (vi,v 2 ), (1) 

where £ is the multipole moment, A c and n are the clustering 
amplitude and index, Cp is the Poisson shot noise, and £q = 
3000 is the pivot scale. This differs from the form adopted 
by Amblard et al. (2011) only in the choice of pivot £q. 
Unlike PI 1 we model the clustering and shot noise as separate, 
independent components rather than modeling their sum as one 
power law. Motivated by the apparent uniformity in angular 
scale dependence across the different spectra (see Figure 1), we 


fit for a single, frequency-independent, value of n. We fit for a 
separate shot-noise level in each of the 16 Planck, BLAST, and 
BLASTxACT spectra (see Section 3.2). 

The SED of the CIB , over the range of frequencies considered, 
has been found to be well described by a modified blackbody of 
the form (e.g., Fixsen et al. 1998; Lagache et al. 1999; Gispert 
et al. 2000) 

/cffi(v) a v^B(v, 7k). (2) 

In reality, the sources making up the CIB lie at a range of 
redshifts, with a variety of luminosities and dust temperatures 
(e.g., Haiman & Knox 2000; Knox et al. 2001; Coppin et al. 
2008; Pascale et al. 2009; Hwang et al. 2010). Equation (2) is 
thus an approximation to the true CIB SED, which consists of 
the sum of many different (approximately) modified blackbody 
spectra, and the quantities ft and T e s are not to be interpreted as 
physical parameters. 

Pll found that the SED of the CIB anisotropies measured 
by Planck is also well described by the Gispert et al. (2000) 
modified blackbody with emissivity spectral index (3 = 1 .4 ± 
0.2 and effective temperature T e g = 13.6 ± 1.5 K. We assume 
that the clustering power SED can also be described in this 
form and parameterize the frequency dependence of the auto- 
spectrum clustering power amplitude as A c = (/ c ) 2 where 

Ic(v) = Io B(y, T eff ), (3) 

with a single emissivity index and effective temperature. 

Different frequencies are sensitive to IR sources at different 
redshifts, with the importance of higher-redshift sources increas- 
ing at lower frequencies (e.g., Haiman & Knox 2000; Knox et al. 
2001; Chapin et al. 2009; Marsden et al. 2009). As a result we 
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would not expect there to be 100% correlation between the 
IR sources in different bands (particularly between the widely 
spaced ACT and BLAST bands) and we therefore include a 
clustering correlation factor / cotr as a free parameter for each 
cross-spectrum. The cross-spectrum clustering amplitude is then 
A c (lh, W) = /corrOl, V2)/ c (Vi)/ c (v 2 ). 

3.1. Shot Noise 


The auto-spectrum shot noise is given by (e.g., Scott & White 
1999) 

f f s “ , d 2 N 

Cp = I dz I dSS 2 -—(S,z), (4) 

J Jo dSdz 

where S cut is the flux cut applied to the map and d 2 N/dSdz are 
the differential source counts. 

In Pll the shot-noise levels were fixed using the IR galaxy 
evolution model of Bethermin et al. (201 1, hereafter B 1 1). This 
model parameterizes the evolution of the galaxy luminosity 
function and is fitted to number counts over a wide range of 
IR wavelengths. While the model broadly fits the available data, 
there are discrepancies. The model underpredicts by ~40% the 
shot-noise levels measured at 500 /x m by BLAST (Viero et al. 
2009) and 220 GHz by SPT (Hall et al. 2010). We do not 
apply any priors on the shot-noise levels when constructing our 
template. We may, however, expect a considerable degeneracy 
between the shape of the clustering component and size of 
the shot noise for the Planck data, because Planck is not 
able to probe the small scales where the shot noise becomes 
dominant. As a result, we consider the effect of adopting the 
B 1 1 model predictions as priors on the Planck shot-noise levels 
in Section 4.1.1. 

In terms of the source flux and S c ut we can write the clustering 
power as (e.g., Tegmark et al. 2002; Viero et al. 2009) 

/ ( d V \ — i / d c \ 2 
dz \dz {z 7 l^ (z) ) P ^ k ’ z ^/xizb (5) 


where / is the comoving distance, dV /dz = y 2 dy/dz is 
the comoving volume element, P ga ] is the IR galaxy power 
spectrum, and the redshift distribution of the flux, dS/dz, is 
given by 


dS r 

7 2) = 1 


Scm d 2 N 
d-S S——{S, z). 
dSdz 


(6) 


The factor of S 2 in Equation (4), compared to idS /dz.) 2 in 
Equation (5), suggests that the removal of the highest-flux 
sources will have much less impact on the clustered power than 
on the shot noise. The B 1 1 model predicts that, for instance, 
applying a flux cut of 250 mJy (the cut applied to the BLAST 
350 /xm map in H12) to the Planck 857 GHz map would result 
in a ~10% reduction in the shot-noise level compared to the 
Planck flux cut of 710 mJy, but that (dS/dz) 2 would be reduced 
by <1% at z ~ 0.2 and virtually unaffected for z > L We 
allow for the dependence of the shot-noise levels on flux cut by 
fitting separate shot-noise levels for each spectrum, as described 
in the next section. We assume that the effect on the clustering 
power from applying different flux cuts is negligible for the data 
considered. Studies using higher-resolution Herschel maps will 
be able to investigate the dependence of the clustering power on 
flux cut in more detail. 


3.2. MCMC Fitting 

We perform a simultaneous fit to the seven auto-spectra (four 
Planck and three BLAST) plus the nine cross-spectra (three 


BLAST x BLAST and six BLAST x ACT). The model spectra 
are binned for comparison to each data spectrum, with likelihood 

16 

- 2ln L(d\9) = [ c ? 3ta - Cj nodel (0)] 2 /o r ? , (7) 

i=i 

for model parameters 9, data vector C‘ klta for the zth spectrum, 
and binned model spectra C™ odel (0). Covariances between the 
16 spectra are neglected. 

We find that the clustering SED parameters / n , /l, and 
T e ff (Equation (3)) are highly degenerate. Only two of these 
parameters are really independent and so we fix T e g- to the best- 
fit value, 9.7 K. We choose Vq = 530 GHz to minimize the 
degeneracy between Iq and /3 . 

We therefore fit for 37 parameters: «, the clustering index, Iq, 
fi, 16 shot-noise levels (one for each spectrum), nine cross- 
spectrum correlation factors (one for each cross-spectrum), 
and nine photometric calibration parameters (one for each 
Planck, BLAST, and ACT band). For each calibration factor 
we enforce a Gaussian prior centered at unity, with spread given 
by the nominal uncertainty listed in Section 2. The covariance 
between the BLAST calibration factors from Truch et al. (2009) 
is included. Apart from the calibration factors, all priors are 
uniform. 

Given the high dimensionality, we estimate the posterior prob- 
ability distribution using Markov Chain Monte Carlo (MCMC) 
analysis (Metropolis et al. 1953), using the sampling methods 
and convergence test described in Dunkley et al. (2005). Chains 
used for analysis are about 10 6 steps in length, and are used to 
calculate one-dimensional marginalized parameter values and 
errors. 

To judge the goodness of fit of the model, we are also 
interested in the maximum likelihood. However, the peak of the 
likelihood distribution occupies only a small part of this high- 
dimensional parameter space, so the minimum y 2 sampled in 
a chain of about a million steps is significantly larger than the 
true value (A/ 2 ~ 10 for a simulated 37-dimensional Gaussian 
distribution). There are numerous statistical methods to find the 
true peak of a distribution. We adopt a simple modification to the 
Metropolis algorithm in which chains start from the best-fitting 
point and then make a step in parameter space only when the 
posterior is improved, using a reduced trial step size. We have 
tested this prescription on simulated Gaussian distributions, and 
the peak is found to within A/ 2 = 0.1 in ~2 x 10 4 steps for 
37 dimensions. We emphasize that this modified code is used 
only to assess the global goodness of fit, not to estimate the 
marginalized parameter distributions. 

To estimate the model spectra at each frequency, it is not 
sufficient to use the nominal values of v for the BLAST and 
ACT bands in Equation (3), due to the finite width of each filter. 
The correct effective v values depend on the SED of the emission 
mechanism. The clustering SED is estimated iteratively by 
repeating the MCMC fitting; each time the best-fit clustering 
SED is integrated through the BLAST and ACT filter profiles 
and an improved estimate of the effective v values obtained. 
We found that after three iterations the SED converged, with 
effective frequencies of 1248, 829, and 607 GHz (effective 
wavelengths 240, 362, and 494 /xm) for the nominal BLAST 
250, 350, and 500 /xm bands, and 220 and 150 GHz for the ACT 
218 and 148 GHz bands. While the clustering SED is rising 
steeply at the ACT bands, the shift from nominal frequency is 
reduced by the narrowness of the ACT filters. For the Planck 
data the nominal frequencies are used, because the temperature- 
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Table 1 

Power-law Clustering Template 


Parameter 

Value 

n a 

1.25 ±0.06 

/oTysr- 1 ) 3 

(2.43 ±0.16) x 10~ 9 

p 

2.20 ± 0.07 

Tat (K) 

9.7 

vo (GHz) 

530 

to 

3000 

Z 2 /dof 

132/122 


Note. a A strong (89%) anti-correlation exists be- 
tween n and Iq. 


to-flux conversion factors provided in Pll are given such that 
the flux is correct at the nominal frequencies. 

4. RESULTS 

The model fits the data well. We find a best- fit / 2 of 132 for 
122 degrees of freedom (150 data points minus 28 parameters; 
we have not counted the calibration parameters since they are 
strongly constrained by priors), giving reduced-/ 2 = 1.08. 
The marginalized mean values of the clustered IR template 
parameters are given in Table 1. No errors are given for T e g, Vq, 
or £o since these parameters are fixed as described in Section 3.2. 
Figure 2 shows the best-fit clustered IR source power for each 
of the bands included in the fit. The best-fit shot-noise levels 
have been subtracted. The y 2 contribution from each individual 
spectrum is also shown; in addition to the total / 2 being good 
we find that there are no spectra which are not individually well 
fitted by the model. 

We plot the BLAST 350 and 500 /xm spectra on the same axes 
as the Planck 857 and 545 GHz spectra. The BLAST spectra 
are color-corrected to account for the difference in the bandpass 
filters by taking the ratio of our clustering power predictions 
for the relevant Planck and BLAST bands. We find these color 
correction factors to be 1.07 and 0.63 for the 350 and 500 /xm 
BLAST spectra, respectively. Pll found values of 1.05 and 0.7 
by integrating the SED of Gispert et al. (2000) through the 
BLAST and Planck filters. These values differ slightly because 
our best-fit clustering SED is slightly different — see Figure 3 
for a comparison. 

The clustering power at l = 3000 is shown as a function 
of frequency in Figure 3. Also shown is the scaled Gispert 
et al. (2000) SED (dashed line), which was fitted to FIRAS and 
DIRBE measurements of the CIB intensity spectrum. It falls 
off more slowly with decreasing frequency than our clustering 
SED. Some difference between the CIB mean and anisotropy 
SED may be expected, since the contribution of a source to the 
anisotropy SED depends on its clustering as well as spectral 
properties. 

IR clustering power predictions at t = 2000 and 3000, 
calculated from our template for several Planck, ACT, and 
SPT bands, are given in Table 2. The conversion from flux 
to temperature units was calculated in each case by integrating 
7 C from Equation (3) through the relevant bandpass filter. The 
effective frequencies (i.e., single frequency value that gives the 
same clustering amplitude as integrating through the filter) are 
also given. Since the IR clustering amplitude is rising strongly 
with frequency at the CMB bands (see Figure 3), the different 
filter profiles lead to significant differences in the clustering 
power even for bands with closely spaced nominal frequencies. 


Table 2 

IR Clustering Predictions from Our Template - 1(1 + l)C£ lust /2jr (/rK 2 ) 


Band 

Veff 

(GHz) 

l = 2000 

i = 3000 

Planck 100 GHz 

104 

0.49 ± 0.07 

0.67 ± 0.09 

Planck 143 GHz 

146 

2.9 ± 0.3 

3.9 ± 0.4 

Planck 111 GHz 

226 

45 ±4 

61 ±6 

ACT 148 GHz 

150 

3.2 ± 0.4 

4.4 ± 0.5 

ACT 218 GHz 

220 

36 ±3 

49 ±4 

SPT 95 GHz 

99 

0.38 ± 0.05 

0.51 ±0.07 

SPT 150 GHz 

156 

4.1 ±0.4 

5.5 ± 0.6 

SPT 220 GHz 

221 

37 ±3 

50 ±5 


The amplitude of the clustered power may vary by a factor of 
ten across a single filter, as shown in Figure 3. 

The predictions at l = 3000 can be compared to the 
ACT results from Dunkley et al. (2011) and SPT results from 
Shirokoff et al. (2011). Dunkley et al. find a best-fit clustering 
amplitude of 4.6 ±1.1 and 54 ± 13 /xK 2 at 148 and 218 GHz, 
respectively (we have calculated these uncertainties by adding 
their statistical and systematic error estimates in quadrature). 
Shirokoff et al. find 6.1 ± 0.8 and 57 ± 8 /xK 2 at 150 and 
220 GHz, respectively, for their base-line model, which includes 
a power-law IR clustering component with index n = 1.2. 
Our predictions are in excellent agreement with these results. 
Figure 4 shows ACT and SPT data at 220 GHz with our 
clustering template predictions overplotted, along with an IR 
shot-noise component and a ACDM CMB power spectrum. 

We also make predictions for the spectral index a of the 
clustered component. At the CMB frequencies the clustering 
SED can be approximated as a power law, I c (v) oc v“ ; however, 
the CMB bands do not lie strictly in the Rayleigh-Jeans limit 
of the modified blackbody adopted in our model. Consequently 
the equivalent power-law slope at 150 GHz differs from that at 
220 GHz: we find q! 15 o = 3.78 ± 0.07 and a 220 = 3.56 ± 0.07. 
We also calculate an effective spectral index between 150 and 
220 GHz, ai 5 o_ 22 o = 3.68 ± 0.07. This result is consistent 
with ACT and SPT findings of 3.69 ± 0.14 and 3.58 ± 0.08, 
respectively. 

We have considered only the frequency dependence of the 
IR clustering component, not the shot noise. Investigating shot- 
noise frequency dependence is more difficult, due to stronger 
flux cut dependence and the fact that the Planck data are limited 
to large scales; future Herschel, Planck, ACT, SPT, and cross- 
correlation data will help overcome these issues. 

4.1. Validating Assumptions 

While our single-index power-law clustering model provides 
a good fit to the data in terms of y 1 , in this section we attempt 
to further validate our assumptions regarding the scale and 
frequency dependence of the clustered power. 

4.1.1. Shot Noise 

Table 3 shows the one-dimensional marginalized values 
of the shot-noise levels from our MCMC chains (Cp from 
Equation (1)). It should be noted that, as expected, the Planck 
shot-noise levels and the clustering index n are highly correlated. 
The shot-noise predictions for each auto-spectrum from the B 1 1 
model are given for comparison. Our values are consistent 
with the model within 1.5er in all cases. We repeated the 
MCMC analysis described in Section 3.2 using the B1 1 model 
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Figure 2. IR point-source clustering power spectra from Planck (diamonds), BLAST (squares), and BLAST/ACT cross-correlations (triangles). The solid lines are 
the best-fit power law with scale dependence C| lust <x l ~ n — we find n = 1.25 ± 0.06. The frequency dependence of the clustering SED is described by a modified 
blackbody (see Equation (3) and Table 1). Poisson shot noise has been subtracted for each panel (shot-noise levels from our fitting are given in Table 3). For the 
combined fit the / 2 /degree of freedom is 132/122. The contribution to the total x 2 from the individual spectra is included in each panel. Color-correction factors of 
1.07 and 0.63 have been applied to the BLAST 350 and 500 /j m spectra so they are directly comparable with the Planck 857 and 545 GHz spectra, respectively, as 
described in the text. 

(A color version of this figure is available in the online journal.) 


predictions as Gaussian priors on the Planck shot-noise levels 
and found a best-fit / 2 /dof = 139/122, an increase of A/ 2 = 7 
compared to the case with uniform priors. There was minimal 
change in the clustering template parameters (e.g., n = 1.24 ± 
0.03 with the shot-noise prior). This suggests that uncertainty in 
the Planck shot-noise levels is not having a significant impact 
on our results. 

4.1.2. Photometric Calibration 


along with the nominal uncertainties. We would expect the 
calibration parameters to be consistent with unity within the 
nominal uncertainty. If this were not the case it could indicate 
that the modified blackbody spectrum in Equation (3) was 
not a suitable description for the frequency dependence of 
the clustering amplitude; however all are consistent with the 
nominal values. 

4.1.3. Planck Beam Uncertainty 


The marginalized values for the photometric calibration 
parameters from our MCMC chains are shown in Table 4 
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As stated in Section 2, we neglected any possible correlation 
in the Planck beam uncertainty across different angular scales. 
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Figure 3. Frequency dependence of IR point-source clustering power at 
i = 3000. We assume this frequency dependence can be described by a 
modified blackbody: Q^Q 00 (y) oc [v@ B(v, T e ff)] 2 . The solid lines show the 
best-fit and 1 a uncertainties from our fitting; parameter values are given in 
Table 1 . The dashed line is the SED of Gispert et al. (2000), which has been 
scaled for comparison. The rectangles show the Planck , BLAST, ACT, and SPT 
bandpass filter FWHM values in the horizontal direction. The vertical extent 
of the rectangles shows the variation of the clustering power across each filter; 
in some cases, this variation is a factor of ten or more due to how steeply the 
clustering SED rises with frequency. The shaded rectangles show the bands 
that were used for the fitting in this paper. Spectra from the other bands are 
either not yet available or were not used, due to the risk of biasing results by 
assumptions regarding the separation of the CMB and other components from 
the IR contribution. 


(A color version of this figure is available in the online journal.) 



Figure 4. Angular power spectra at 220 GHz measured by ACT (Das et al. 
2011) and SPT (Shirokoff et al. 2011), with a theoretical model for CMB and IR 
point sources overplotted. The clustered IR component is given by the power- 
law model fit to Planck , BLAST, and BLAST x ACT data in this analysis (see 
Equations (1) and (3) and Table 1 for parameter values). The upper and lower 
dashed lines correspond to la error bounds. The lensed CMB spectrum is that 
of the ACDM model with parameters derived from WMAP (Komatsu et al. 
2011). An IR shot-noise component of size t{t + l)Cp/27r|£ = 3ooo = 78 pK 2 
(consistent with ACT and SPT measurements) has also been plotted. Additional 
radio source and kinetic SZ power are subdominant at this frequency and are 
not included. 

(A color version of this figure is available in the online journal.) 


We re-ran the MCMC chain enforcing a 100% correlation 
in the beam uncertainty in different €-bins for each Planck 
spectrum (introducing off-diagonal elements in the likelihood — 
Equation (7)) and found minimal change in the marginalized 
parameter values (<0.3cr in all cases). While the Planck beam 


Table 3 

One-dimensional Marginalized Shot-noise Levels 


Band 

Marginalized Value 

Bethermin et al. (201 1) 


(Jy 2 sr” 1 ) 

(Jy 2 sr- 1 ) 

Planck 857 GHz 

2200 ± 2500 

5920 ± 370 

Planck 545 GHz 

1700 ± 700 

1 150 ± 90 

Planck 353 GHz 

210 ±60 

138 ± 22 

Planck 111 GHz 

16 ± 6 

12 ±3 

BLAST 250 pm 

11600 ± 2300 

1 1600 ± 2100 

BLAST 350 /zm 

5200 ± 1200 

5050 ± 1080 

BLAST 500 /zm 

1300 ±390 

1680 ± 480 

250 pm x 350 pm 

5700 ± 1400 


250 pm x 500 pm 

1830 ±690 


350 pm x 500 pm 

820 ± 490 


250 pm x 218 GHz 

240 ± 130 


250 pm x 148 GHz 

90 ±70 


350 |im x 218 GHz 

240 ± 70 


350 pm x 148 GHz 

110 ±40 


500 (jm x 218 GHz 

66 ±30 


500 fim x 148 GHz 

55 ± 20 



Table 4 


One-dimensional Marginalized Calibration Parameters 

Band 

Marginalized Value 

Nominal Uncertainty 

Planck 857 GHz 

1.05 ±0.06 

7% 

Planck 545 GHz 

0.99 ± 0.04 

7% 

Planck 353 GHz 

1.01 ±0.02 

2% 

Planck 217 GHz 

1.00 ±0.02 

2% 

BLAST 250 pm 

1.05 ±0.08 

9.5% 

BLAST 350 pm 

1.03 ±0.07 

8.7% 

BLAST 500 /zm 

1.02 ±0.07 

9.2% 

ACT 218 GHz 

1.03 ±0.07 

7% 

ACT 148 GHz 

1.00 ±0.02 

2% 


uncertainties at l ~ 2000 are comparable to or larger than the 
statistical uncertainties, the template parameters are primarily 
constrained by the large-scale Planck data (where the beam 
uncertainties are sub-dominant) and the BLAST data. 


4.1.4. Comparison with Broken Power Law 


We assumed that the scale dependence of the clustering power 
can be described using a single index n. In this section, we 
consider fitting the scale dependence with a broken power law 
of the form 


cf Jst oc 


r n 1 
£-«2 


if l < 4 
if l > lb, 


( 8 ) 


where n\, n 2 , and lb are free parameters. When we repeat the fit 
to the Planck, BLAST, and BLAST x ACT spectra with this new 
scale dependence we find n\ = 1.24 ± 0.06, n 2 = L2 ± 0.2, 
and l h = 1100 ± 300 with a best x 2 /dof = 131/120. This is 
an improvement of only A/ 2 = 0.9 with two additional fitted 
parameters. For the fit with priors on the Planck shot-noise levels 
from the B1 1 model the improvement is only A/ 2 = 0.4. We 
conclude that the data do not show a significant preference for 
a broken power law. 

Keisler et al. (2011) adopted a template of the form given 
in Equation (8) with n\ = 2.0, n 2 = 1.2, and lb = 1500 to 
model the IR clustered power when extracting cosmological 
information from the small-scale SPT 150 GHz CMB power 
spectrum. We find that this form is not a good fit to the Planck 
data, with / 2 /dof = 22/7 when we fit to the Planck 217 GHz 
spectrum with n\, n 2 , and lb fixed to the above values and using 
uniform priors on the clustering amplitude and shot noise. 
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Table 5 

Degree of Clustering Cross-correlation, / corr 


Band 1 

Band 2 

Av 

(GHz) 

Marginalized Value 

250 fim 

350 fim 

340 

0.98 ± 0.30 

250 fim 

500 fim 

600 

1.23 ±0.31 

350 pm 

500 fim 

260 

1.41 ±0.28 

250 fim 

218 GHz 

982 

0.66 ± 1.09 

250 fim 

148 GHz 

1052 

0.30 ± 1.83 

350 fim 

218 GHz 

642 

0.64 ± 0.56 

350 /am 

148 GHz 

712 

0.39 ± 1.05 

500 /am 

218 GHz 

382 

1.82 ±0.47 

500 fim 

148 GHz 

452 

0.79 ± 0.87 


4.1.5. Comparison with Linear Bias Model 

H12 found that the BLAST and BLAST x ACT data are well 
fitted by assuming the IR galaxy power spectrum in Equation (5) 
is given by 

Pgai(k, z ) = b 2 P DM (k, z ), (9) 

where b is the linear bias factor and Pdm is the linear dark matter 
power spectrum, and using the B 1 1 predictions for the redshift 
distribution of the flux, dS/dz. PI 1 also found that the Planck 
data from each band are well fitted by this model (also using the 
B 1 1 dS/dz ) if no prior is enforced on the shot-noise levels. The 
bias levels are not consistent; when fitting for a single, redshift- 
independent value of b, H12 found b — 5.0 ± 0.4 whereas PI 1 
found b — 2. 18 ± 0. 1 1 for the Planck 545 GHz spectrum. 

The linear bias model is strongly rejected by the combined 
Planck and BLAST data even without any shot-noise prior. 
Fitting a linear bias model with single-value bias to the Planck 
857 GHz and BLAST 350 /im spectra together, using the B 1 1 
dS/dz and multiplying the BLAST data points by a color- 
correction factor of 1.07 (see Section 4), yields/ 2 /dof = 39/16. 
For a joint fit to the Planck 545 GHz and BLAST 500 /rm spectra 
we find / 2 /dof = 43/16. This result is driven by the shape 
of the linear matter power spectrum rather than the choice of 
clS/dz : we repeated the fitting using the predictions of Marsden 
et al. (2011) rather than Bll and found y 2 /dof of 36/16 and 
42/16 for the 857 GHz/350 /tm and 545 GHz/500 pm spectra, 
respectively. Neither the Planck nor BLAST data alone were 
able to rule out the linear bias model without a shot-noise prior 
because of the limited angular scales probed. Shirokoff et al. 
(2011) similarly found that the linear bias model could not be 
ruled out using only SPT data from l > 2000. 

4.2. Cross-correlations 

Table 5 shows the marginalized degree of clustering correla- 
tion, /con- for each cross-spectrum included in our fitting, along 
with the separation Av between the two bands. The BLAST x 
BLAST cross-spectra are consistent with 100% correlation. Our 
conclusions regarding the degrees of correlation between the 
ACT and BLAST bands are limited by the data quality. A de- 
crease in correlation with increasing band separation would be 
consistent with the sources lying at a range of redshifts, with the 
higher-redshift sources being of greater relative importance at 
the longer wavelengths (e.g., Haiman & Knox 2000); however 
the current data are not of sufficient quality to confirm this. 

Two of the marginalized mean / corr values lie more than 1 a 
above unity, while none lie more than la below. Measuring 
/con- > 1 may indicate that the angular scale dependence of the 
cross-spectra clustering power is not described by the same 


single-index power law as the auto-spectra clustering, since 
there is no physical explanation for a correlation in excess of 
100%. This could be the case even with no worsening of the / 2 , 
due to the limited angular scales probed by the BLAST x ACT 
data. To test that this is not having a significant effect on our 
clustering template, we repeated the MCMC fitting described in 
Section 3.2 using only the auto-spectrum data. We found that 
the values of n, [i, and Iq change by <0.5a compared to the 
fit with the cross-spectra included. We conclude that while the 
data may be hinting that the cross-spectrum clustering power 
has a different shape to the auto-spectrum clustering, this is not 
significantly biasing our template. 

Further submillimeter-millimeter cross-correlation studies 
(e.g., Herschel x ACT/SPT) are clearly required to provide 
more insight into the distribution of the sources with redshift, 
to constrain the angular scale dependence of the cross-spectrum 
clustering, and to investigate how the clustering shape changes 
with increasing band separation. 

5. CONCLUSIONS 

We have found that a power-law model for the IR point-source 
clustering is adequate to simultaneously fit Planck, BLAST, 
and cross-correlated BLAST/ACT power spectrum data over 
a broad range of frequency (150 GHz < v < 1200 GHz) 
and angular scale (multipole moment 100 < l < 9000). We 
find that the clustering power varies with angular scale as l~ n 
with n — 1.25 ± 0.06 and that the SED of the clustering can 
be described as a modified blackbody with emissivity index 
/ = 2.20 ± 0.07 and effective temperature 7/. rf - = 9.7 K. 
Our work does not rely on any assumptions regarding the 
physical properties of the IR sources (host halos, redshift 
distribution, etc.). 

As well as providing a simple template for use in CMB 
foreground subtraction, we have established that the Planck 
and BLAST/BLAST x ACT data sets appear compatible when 
bandpass filters, flux cut, and calibration are accounted for, as 
described in Sections 2 and 3. We make predictions for the IR 
clustering power for the Planck, ACT, and SPT CMB bands; 
our predictions for the ACT and SPT bands at around 150 
and 220 GHz are fully consistent with existing measurements 
(Dunkley et al. 2011; Shirokoff et al. 201 1). 

There is uncertainty in the Planck shot-noise levels because 
Planck does not probe small enough scales for the shot noise 
to become dominant. We find reasonable consistency between 
the shot-noise levels from our fitting and the predictions of the 
parametric IR galaxy evolution model of Bethermin et al. (20 11). 
We repeated our fitting using this model’s predictions as priors 
on the Planck shot-noise levels and found that there was minimal 
effect on the scale or frequency dependence of the clustering 
power. Number counts obtained via P(D) analysis (probability 
of deflection — modeling the probability distribution function of 
observed flux in each map pixel — e.g., Patanchon et al. 2009; 
Glenn et al. 2010) may be useful for constraining shot-noise 
levels in future work. 

Upcoming data from Herschel, Planck, ACT, SPT, and 
cross-correlations will allow us to look for any variations 
in the clustering angular scale dependence with frequency. 
Understanding this scale dependence in the millimeter CMB 
bands is important not only for extracting unbiased estimates 
of cosmological parameters, but also for constraining other 
components of the measured spectrum, such as the thermal 
and kinetic SZ effect, and any cross-component correlations, 
for instance SZ-IR. 
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Over the range of angular scales considered we would 
expect to be probing both the linear and nonlinear clustering 
regimes. This is supported by the fact that the combined 
Planck and BLAST data strongly reject a linear bias clustering 
model (Section 4.1.5). The linear and nonlinear components 
have different scale dependence, with the linear component 
dominating on large scales and the nonlinear on small scales. 
The fact that our study has found that a power-law scale 
dependence provides a good fit to the Planck and BLAST 
IR point-source clustering power is not inconsistent with this 
picture; instead it suggests that the sum of the linear and 
nonlinear components is sufficiently close to a power law that 
the current data are unable to reveal any deviations. We can 
draw a comparison to measurements of the angular correlation 
function, w{9), of luminous red galaxies and galaxies in the main 
SDSS sample, which have revealed deviations from power-law 
behavior (e.g., Zehavi et al. 2004, 2005; Blake et al. 2008), 
whereas earlier studies of galaxy clustering were unable to do 
so (e.g., Zehavi et al. 2002). We may expect similar deviations 
to be detected in future IR source power spectra, in particular 
with Herschel /SPIRE data, which will yield higher quality data 
than BLAST out to smaller scales with its improved angular 
resolution. 

The correlation function of resolved submillimeter sources 
(e.g., Cooray et al. 2010; Maddox et al. 2010) provides comple- 
mentary information to the power spectra of unresolved sources 
considered in this work. Current data are limited; Guo et al. 
(2011) measure £(r) oc r~ Y with y ~ 2, corresponding to 
C £ clust oc £“ 1 , for low-redshift (z < 0.5) Herschel - ATLAS galax- 
ies, although their result is consistent with ours within errors. 
The full Herschel - ATLAS survey (Eales et al. 2010) will cover 
30 times more sky and provide far tighter constraints. Discrep- 
ancies between the C£ lust and f (r) or w(Q) measurements would 
indicate differences in clustering properties of the bright sub- 
millimeter galaxies compared to the fainter population, and thus 
both statistics are important for constraining models of galaxy 
clustering and evolution. 
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